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SECTION  I 
INTRODUCTION 

HULL  Is  a computer  code  that  solves,  within  an  Eulerian  mesh,  the  hydro- 
dynamic  equations  of  continuity,  momentum,  and  energy  with  a finite  difference 
scheme  (ref.  1).  The  code  has  an  equilibrium  radiation  diffusion  treatment. 
However,  the  assumption  of  equilibrium  diffusion  severely  limits  the  types  of 
problems  that  can  be  addressed.  Moreover,  the  manner  in  which  the  equilibrium 
diffusion  is  implemented  requires  an  unacceptably  small  timestep.  For  these 
reasons  we  have  incorporated  a treatment  of  non-equilibrium  radiation  diffusion 
that  includes  a flux-limiter  for  optically  thin  regions. 

Non-equilibrium  diffusion  relaxes  the  assumption  that  the  material  and 
radiation  must  have  the  same  temperature.  However,  we  have  not  included  a multi - 
frequency  treatment,  so  we  are  still  forced  to  assume  the  radiation  field  has  a 
black  body  distribution,  although  the  black  body  temperature  is  not  required  to 
be  equal  to  the  material  temperature.  The  flux-limiters  have  been  added  to 
extend  the  diffusion  theory  to  optically  thin  regions.  The  diffusion  coeffi- 
cient is  modified  so  that  the  flux  goes  over  to  the  free-streaming  results  when 
the  mean-free-path  becomes  long. 
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SECTION  II 

THE  DIFFUSION  EQUATION  IN  CYLINDRICAL  COORDINATES 
We  write  the  diffusion  equation  in  cylindrical  coordinates  as 


and  the  associated  change  in  the  material  energy  density  as 

3t p,cac(aT  * E) 

Here  E is  the  radiation  energy  density,  p is  the  material  density,  <_  is  the 

a 

absorption  opacity,  c is  light  velocity,  a is  the  radiation  constant,  T Is  the 
material  temperature,  Em  Is  the  material  energy  density,  and  FR  and  F^  are  the 
radiation  fluxes  In  the  R and  Z directions,  respectively.  This  equation 
incorporates  the  non-equilibrium  diffusion  approximation;  the  radiation  energy 
density  is  not  assu"  %o  be  aT^,  the  black  body  energy  density.  We  note  that 
no  material  radiation  energy  exchange  due  to  Compton  scattering  is  included. 
While  we  do  not  incorporate  the  Compton  energy  exchange,  Compton  scattering  is 
Included  In  the  diffusion  coefficients  and  hence  In  the  computation  of  the 
radiation  fluxes. 

The  material-radiation  energy  exchange  Is  done  with  the  operator  splitting 
method  as  discussed,  for  example,  by  Richtmyer  (ref. 2).  Operator  splitting  Is 
also  used  to  solve  the  diffusion  equation;  we  do  the  R diffusion  separately 
from  the  Z diffusion. 

FR  and  F^  are  taken  to  be 


F 


R 


f 


cX  3E 
R T 


2.  Richtmyer,  R. , and  Morton,  K.,  Difference  Methods  for  Initial-Value 
Probl ems , Interscience  Publishers,  New  York,  1967. 
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c — f CX  3E 
FZ"fZl9Z 


where  X is  the  total  mean  free  path  (absorption  and  scattering),  and  f^  are  fj 
are  the  "flux-limiters"  in  the  R and  Z directions,  respectively,  f^  and  f^ 
are  computed  as 


1 + 3 X 


3 Jin  E 
3R 


(l  + 3 exp  | - j X 


3 Jin  E 

3R 


1 I + 3 exp  j-  ^-X  | 


3 ^n  E 


We  see  that  our  expressions  for  and  are  the  standard  diffusion  fluxes 
except  for  the  flux  limiters.  For  small  X (when  diffusion  should  apply),  the 
flux  limiters  are  near  unity,  and  we  have  ordinary  diffusion  theory.  When  X 
becomes  sufficiently  large  (which  depends  on  the  scale  of  the  problem)  so  that 
the  diffusion  approximation  breaks  down,  the  flux  limiters  reduce  the  diffusion 
coefficients  so  that  the  radiation  flux  cannot  exceed  the  product  of  light 
velocity  and  the  energy  density,  which  is  the  maximum  physically  allowable  flux. 
To  see  this  limiter  behavior,  consider  F^  with  X large.  We  neglect  the 
exponential  term  since  it  has  a negative  argument  of  large  absolute  value.  Then 
we  have 


- cX  / 

3 ll  + 


C - CX 3E 

R , . . 1 , 3E  I 3R 

3 + x r I 3R  I 

Assuming  X is  sufficiently  large  so  that  we  may  neglect  the  3 in  the  denominator, 
we  have 
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■ = 3E/9R 

R ' " 7aF/3R| 


cE 


That  is,  | Fp | * cE,  with  the  direction  of  flow  determined  by  the  direction  of 
the  gradient  in  E. 

Our  discussion  has  been  slightly  misleading  in  that  we  have  concentrated  on 
the  value  of  X.  Of  more  importance  is  the  size  of  X compared  with  where 


tE  is  a radiation  scale  length  (or  scale  height).  If  X » *E,  then  the  flux 
limiters  become  operational.  However,  if  » X,  fR»  fz  * 1,  and  we  have 
normal  diffusion. 

A more  complete  description  of  flux  limited  diffusion  theory  is  given  by 
Winslow  (refs.  3 and  4).  As  indicated  by  Winslow,  much  of  the  flux  limiter 
development  has  been  done  by  LeBlanc  and  Wilson  of  the  Lawrence  Livermore 
Laboratory. 


3.  Winslow,  A.M. , Improved  Flux  Limiter  for  Asymptotic  Neutron  Diffusion  Calcu- 
lations, UCIR-378,  Lawrence  Livermore  Laboratory,  Livermore,  CA,  April  1969. 

4.  Winslow,  A.M. , "Extensions  of  Asymptotic  Neutron  Diffusion  Theory," 

Nuclear  Science  and  Engineering,  32,  pp  101-110,  1968. 
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SECTION  III 
DIFFERENCE  EQUATIONS 

We  adopt  operator  splitting  methods  (e.g.,  ref.  2)  to  separate  our  diffu- 
sion equation  into  a radial  component  and  an  axial  component,  each  of  which  is 
separately  solved  numerically.  Half  of  the  radiation-material  coupling  is  done 
with  each  radiation  sweep.  That  is,  we  assume 


9E  _ 3E  | 3E  | 

at  " at  'R  at  'z 


3Em  aE  aE 

m _ m i , m , 

at  " at  'r  at  'z 


where 


at  (R  = ' R |r  (RFR}  + ?P<ac(aT4  - E) 


H >Z  = " fz  FZ  +!p<ac(al4  ' E> 


aEm 

— I 
at  'R 


at  lz  = ■ 2 pKac ( aT^  - E) 


We  first  derive  the  radial  difference  equation  in  detail.  The  fully 
implicit  difference  equations  are 


E - E« 


- At  /f 
' rar  \ 


E+  " E _ Rex. 

+1/2  R+  " R 3 -1/2 


E - E \ pic  cAt  . 

nr)+  -h-  <aT  - E> 


pK  CAt  . 

E™  - * - -V-  <aT  - E> 
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where  the  E°,  E°,  and  TQ  denote  the  quantities  at  the  beginning  of  the  current 
time  step.  Unless  otherwise  indicated,  the  quantities  are  cell  centered.  Here 
+ (-)  refers  to  the  adjacent  cell  center  at  larger  (smaller)  radius.  We  use 
+1/2  (-1/2)  to  indicate  the  boundary  between  the  current  cell  and  the  + (-) 
cell.  Note  that  AR(no  subscript)  is  computed  as  AR  = R+1/2  " R-l/2’  These 
equations,  coupled  with  the  equation  of  state,  completely  determine  E,  Em,  and 
T.  In  order  to  calculate  these  quantities,  we  first  make  the  approximation 

Em  - e;  ■ ecv<T  - T0> 

where  Cy  is  the  material  specific  he?t.  We  are  then  faced  with  solving  a 4th 
order  system  in  T.  We  choose  to  solve  this  system  by  linearizing  the  4th 
order  term. 

The  linearization  is  taken  to  be 

T4  = T4  + 4T?  (T  - T ) 

0 o'  o' 

= 4T^T  - 3T4 

O 0 

Then 

E • E°  * ■ rarlrn  I+j/2(e+  • E) 


and 

PCV(T  - Tq)  = 
Collecting  terms,  we  have 


( 

i 


6 


Now  define  6 as 


or  solving  for  E 


Concentrating  on  the  material  interaction  terms,  we  have 
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SECTION  IV 


TEST  PROBLEMS 

We  have  run  two  classes  of  test  problems.  First,  we  have  calculated  the 
propagation  of  plane  waves  parallel  to  the  z axis  and  compared  the  results 
with  an  established  one-dimensional  code.  Second,  we  have  placed  a spherical 
source  in  the  center  of  the  two-dimensional  mesh  and  followed  the  radiation 
front  to  test  the  symmetry  characteristics  of  the  difference  equations. 

As  an  example  of  the  plane  wave  problems,  we  present  results  from  a calcu- 
lation with  40  cells  in  the  z direction.  The  source  region  consisted  of  the 

first  five  rows  with  energies  of  2 x 10^4  ergs/gm  while  the  remainder  of  the 

12 

mesh  was  filled  with  an  ambient  energy  of  10  ergs/gm.  Reflective  boundary 

conditions  were  used  on  all  boundaries.  Cells  were  250  cm  square,  and  the 

mean  free  path  was  7500  cm  (30  cells).  Circles  with  dots  in  figure  1 indicate 
the  flux  for  each  of  the  first  28  zones  after  20  cycles.  For  our  comparison 
we  have  chosen  a one-dimensional  code  that  employs  a similar  radiation  flux- 
limited  diffusion  scheme  (ref.  5).  Identical  initial  conditions  were  used  in 
this  code,  and  the  results  are  shown  as  x's  in  figure  1.  The  agreement  is 
excellent  - the  slight  differences  in  the  flux  for  the  outer  cells  is  the 
result  of  using  a different  procedure  to  calculate  an  average  opacity  for  the 
outer  boundary  in  the  two  codes. 

The  second  case  used  a mesh  of  100  x 100  cells  with  cell  dimensions  equal 
to  50  cm.  Our  spherical  source  with  an  energy  density  of  2 x 10^  ergs/gm 
had  a radius  of  25  zones.  The  mean  free  path  was  fixed  at  1000  cm  (20  cells). 
Figures  2 and  3 show  the  initial  conditions  of  this  calculation  for  material 
energy  density  and  radiation  energy  density.  Figures  4,  5,  6,  and  7 illustrate 
the  expanding  source. 


5.  Alme,  M.L.  and  Wilson,  J.R.,  "Numerical  Study  of  Accretion  onto  a Neutron 


Star,"  Astrophysical  Journal , 186,  p 1015,  1973. 
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Figure  1.  Comparison  of  HULL  Non-Equilibrium  to  1-0  Code 


10 


i 


ALTITUDE  M 


10  20 


ENERGY 
40  50 


COOL 
E = 1012 


WARM 
E = 2- 10 ' 


CONTOUR  SCALE 
ERGS/GM 

2 2.500E  + 13 

3 5.000E+13 

4 7.500E+13 

5 1.000E+14 

6 1.25QE+14 

7 1 . 500E  4 14 

8 1.750E  + 14 

9 2.000E  +14 

DXI  * 5.000E  +01 
MIN  = 1 ,OOQE  +12 
MAX  = 2 .OOGE  +14 


.0  5.0  10.0  15.0  20.0  25.0  30.0  35.0  40.0  45.0  50.0 

RADIUS  M 

AFWL  NONEQUILIBRIUM  RADIATION  CALCULATION 
TIME  0.000  NSEC  CYCLE  0.  PROBLEM  23.0071 


Figure  2.  Initial  Configuration  for  Material  Energy  Density 
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Figure  3.  Initial  Configuration  for  Radiation  Energy  Density 
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Figure  4 Material  Energy  After  150  Cycles 
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Figure  5.  Radiation  Energy  Density  After  150  Cycles 
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Figure  6.  Material  Energy  Density  After  250  Cycles 
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Figure  7.  Radiation  Energy  Density  After  250  Cycles 
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